consolidate all rad-to-miles calculations, add 'points' option to arc filter

diff --git a/gpsbabel/arcdist.c b/gpsbabel/arcdist.c
index 7ee828099d12db2fb6b7f5be9302f99a846b3459..1f5f1d2f1ffaf9de3a41eb0b7f7ae6dbb729e302 100644 (file)
--- a/gpsbabel/arcdist.c
+++ b/gpsbabel/arcdist.c
@@ -30,6 +30,7 @@ static double pos_dist;
 static char *distopt = NULL;
 static char *arcfileopt = NULL;
 static char *exclopt = NULL;
+static char *ptsopt = NULL;
 
 typedef struct {
        double distance;
@@ -42,6 +43,8 @@ arglist_t arcdist_args[] = {
        {"distance", &distopt, "Maximum distance from arc", 
                ARGTYPE_FLOAT | ARGTYPE_REQUIRED},
        {"exclude", &exclopt, "Exclude points close to the arc", ARGTYPE_BOOL},
+       {"points", &ptsopt, "Use distance from vertices not lines", 
+               ARGTYPE_BOOL},
        {0, 0, 0, 0}
 };
 
@@ -78,17 +81,24 @@ arcdist_process(void)
            if ( argsfound != 2 && strspn(line, " \t\n") < strlen(line)) {
                 warning(MYNAME ": Warning: Arc file contains unusable vertex on line %d.\n", fileline );
            } 
-           else if ( lat1 != BADVAL && lon1 != BADVAL &&
-                lat2 != BADVAL && lon2 != BADVAL ) {
+           else if ( lat2 != BADVAL && lon2 != BADVAL &&
+                (ptsopt || (lat1 != BADVAL && lon1 != BADVAL ))) {
              QUEUE_FOR_EACH(&waypt_head, elem, tmp) {
 
                waypointp = (waypoint *)elem;
-               dist = linedist(lat1, lon1, lat2, lon2, 
+               if ( ptsopt ) {
+                  dist = gcdist( lat2*M_PI/180.0, lon2*M_PI/180.0, 
+                                  waypointp->latitude*M_PI/180.0, 
+                                  waypointp->longitude*M_PI/180.0 );
+               }
+               else {
+                  dist = linedist(lat1, lon1, lat2, lon2, 
                                waypointp->latitude,
                                waypointp->longitude );
+               }
 
                /* convert radians to float point statute miles */
-               dist = (((dist * 180.0 * 60.0) / M_PI) * 1.1516);
+               dist = tomiles(dist);
 
                if ( waypointp->extra_data ) {
                        ed = (extra_data *) waypointp->extra_data;

diff --git a/gpsbabel/grtcirc.c b/gpsbabel/grtcirc.c
index d059e3c84a5eea5bb32eeb321dd9ecd23106b8af..7a0b01c99181b3bc4390ed78dd27f66247c88399 100644 (file)
--- a/gpsbabel/grtcirc.c
+++ b/gpsbabel/grtcirc.c
@@ -34,6 +34,22 @@ static double dotproduct( double x1, double y1, double z1,
   return (x1*x2+y1*y2+z1*z2);
 }
 
+/*
+ * Note: this conversion to miles uses the WGS84 value for the radius of
+ * the earth at the equator.
+ * (radius in meters)*(100cm/m) -> (radius in cm)
+ * (radius in cm) / (2.54 cm/in) -> (radius in in)
+ * (radius in in) / (12 in/ft) -> (radius in ft)
+ * (radius in ft) / (5280 ft/mi) -> (radius in mi)
+ * If the compiler is half-decent, it'll do all the math for us at compile
+ * time, so why not leave the expression human-readable?
+ */
+
+double tomiles( double rads ) {
+    const double radmiles = 6378137.0*100.0/2.54/12.0/5280.0;
+    return (rads*radmiles);
+}
+
 double gcdist( double lat1, double lon1, double lat2, double lon2 ) {
   return acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2));
 }

diff --git a/gpsbabel/grtcirc.h b/gpsbabel/grtcirc.h
index 94c354a08227138b77172755f10dcec8a3d56790..de84aea4c60983ea37f4a0f1a32b66b8b5a18601 100644 (file)
--- a/gpsbabel/grtcirc.h
+++ b/gpsbabel/grtcirc.h
@@ -23,3 +23,5 @@ double gcdist( double lat1, double lon1, double lat2, double lon2 );
 double linedist(double lat1, double lon1,
                double lat2, double lon2,
                double lat3, double lon3 );
+
+double tomiles( double rads );

diff --git a/gpsbabel/position.c b/gpsbabel/position.c
index 46bbfc0f46320991a05fa55912178146dcc13e8a..3aa8b5feface7af367c9aed6eec3542b917adfa0 100644 (file)
--- a/gpsbabel/position.c
+++ b/gpsbabel/position.c
@@ -150,7 +150,7 @@ position_runqueue(queue *q, int nelems, int qtype)
                                   comp[i]->longitude);
 
                /* convert radians to integer feet */
-               dist = (int)((((dist * 180.0 * 60.0) / M_PI) * 1.1516) * 5280.0);
+               dist = (int)(5280*tomiles(dist));
                
                if (dist <= pos_dist) {
                        switch (qtype) {
@@ -273,7 +273,7 @@ radius_process(void)
                                   home_pos->longitude);
 
                /* convert radians to float point statute miles */
-               dist = (((dist * 180.0 * 60.0) / M_PI) * 1.1516);
+               dist = tomiles(dist);
 
                if ((dist >= pos_dist) == (exclopt == NULL)) {
                        waypt_del(waypointp);